###############################################################################
## AUTHOR: ALAN YAN
## DATE LAST UPDATED: 09/27/2020
## PURPOSE: HET EFFECTS VIA INTERACTION (BENCHMARKING WITH MCCONNELL ET AL)
###############################################################################
rm(list = ls())

#### NOTES ####

#(1) create a tidier function for random forest output

#### PACKAGES ####
library(pacman)
p_load(tidyverse,
       estimatr,
       broom,
       haven,
       hrbrthemes,
       cjoint,
       formula.tools, 
       DeclareDesign,
       emmeans,
       lmtest,
       cregg)

#### FUNCTIONS ####

# ggplot theme
theme_shom_alt <- function (base_size = 11, waffle = FALSE) 
{
  ret <- theme_minimal(base_size = base_size, base_family = "Roboto Condensed") + 
    theme(plot.background = element_rect(fill = "#f5f5f2", 
                                         color = NA), 
          panel.background = element_rect(fill = "#f5f5f2",color = NA), 
          legend.background = element_rect(fill = "#f5f5f2",color = NA))
  if (waffle) {
    ret + theme(axis.text = element_blank())
  }
  else {
    ret
  }
}

# lm robust function that includes omitted categories for plotting
lm_robust_omitted <- function(formula, data, clusters) {
  require(formula.tools)
  require(tidyverse)
  require(haven)
  results <- lm_robust(formula = formula,
                       data = data,
                       clusters = clusters) %>%
    tidy() 
  omitted_cat <- model.frame(formula = formula, data = data) %>% 
    select(2:length(.)) %>% 
    sapply(., levels) %>% 
    lapply(., "[[", 1) %>% 
    unlist() %>%
    as.data.frame()
  omitted_cat_df <- cbind(paste0(rownames(omitted_cat), omitted_cat[,1]), 
                          0, 0, 0, 0, 0, 0, 0, 
                          lhs.vars(formula)) %>%
    as.data.frame() %>%
    zap_labels() %>%
    setNames(., names(results))
  results_with_omitted <- rbind(results, 
                                omitted_cat_df)
  results_with_omitted[,c(2:8)] <- sapply(results_with_omitted[,c(2:8)],as.numeric)
  return(results_with_omitted)
}

#### LOAD DATA ####
df <- read_rds("01-yougov-conjoint/01-data/cleaned/yougov-conjoint-cleaned.rds")

#### CLEAN DATA ####
df %>%
  mutate(
    copartisan =  case_when(
      r_party %in% c(1:3) & cje_political_donations == "Primarily donates to Democrats" ~ 1,
      r_party %in% c(5:7) & cje_political_donations == "Primarily donates to Republicans" ~ 1,
      r_party %in% c(1:3) & cje_political_donations == "Primarily donates to Republicans" ~ 0,
      r_party %in% c(5:7) & cje_political_donations == "Primarily donates to Democrats" ~ 0,
      r_party %in% c(1:3) & cje_political_donations == "Donates to both Democrats and Republicans" ~ 0,
      r_party %in% c(5:7) & cje_political_donations == "Donates to both Democrats and Republicans" ~ 0,
    ),
    outpartisan =  case_when(
      r_party %in% c(1:3) & cje_political_donations == "Primarily donates to Republicans" ~ 1,
      r_party %in% c(5:7) & cje_political_donations == "Primarily donates to Democrats" ~ 1,
      r_party %in% c(1:3) & cje_political_donations == "Primarily donates to Democrats" ~ 0,
      r_party %in% c(5:7) & cje_political_donations == "Primarily donates to Republicans" ~ 0,
      r_party %in% c(1:3) & cje_political_donations == "Donates to both Democrats and Republicans" ~ 0,
      r_party %in% c(5:7) & cje_political_donations == "Donates to both Democrats and Republicans" ~ 0,
    ),
    work_binary = (work_binary - mean(work_binary, na.rm = TRUE))/sd(work_binary, na.rm = TRUE),
    power_binary = (power_binary - mean(power_binary, na.rm = TRUE))/sd(power_binary, na.rm = TRUE),
    responsibility_binary = (responsibility_binary - mean(responsibility_binary, na.rm = TRUE))/sd(responsibility_binary, na.rm = TRUE),
    complaints_binary = (complaints_binary - mean(complaints_binary, na.rm = TRUE))/sd(complaints_binary, na.rm = TRUE),
    r_party = r_party - mean(r_party, na.rm = TRUE)
  ) -> df

#### ANALYZE DATA (NO CONTROLS) - all data ####
#Do people prefer firms with these features?
main_work_td <- lm_robust_omitted(work_binary ~ 
                                    copartisan +
                                    outpartisan +
                                    as.factor(cje_corporate_gov) +
                                    as.factor(cje_corporate_resp) +
                                    as.factor(cje_firm_size) +
                                    as.factor(cje_gender_ownership) + 
                                    as.factor(cje_healthcare) + 
                                    as.factor(cje_hours) + 
                                    as.factor(cje_training) + 
                                    as.factor(cje_location) + 
                                    as.factor(cje_sick_leave) + 
                                    as.factor(cje_parental_leave) + 
                                    as.factor(cje_race_ownership) + 
                                    as.factor(cje_retirement) + 
                                    as.factor(cje_income) +
                                    as.factor(cje_type_of_work) +
                                    as.factor(cje_union) +
                                    as.factor(cje_work_from_home) +
                                    as.factor(cje_team_work) +
                                    as.factor(cje_work_culture)
                                    ,
                                  data = df,
                                  clusters = df$rid
) %>%
  filter(term != "(Intercept)") %>%
  mutate(term_clean = factor(gsub(x = term,pattern = "as.factor",replacement = "")),
         estimate = ifelse(is.na(estimate), 0, estimate),
         std.error = ifelse(is.na(std.error), 0, std.error),
         conf.low.90 = estimate - 1.64*std.error,
         conf.high.90 = estimate + 1.64*std.error) 
main_work_td

#our estimates
esop.estimate <- 0.2016158916 # ESOPs
codeterm.estimate <- 0.1486489704 # codetermination
election.estimate <- 0.0809887104 # workers elect managers
copartisan.estimate <- 0.1065834301 # copartisan
outpartisan.estimate <- -0.2549901041 #outpartisan

#mcconnell et al estimates
mcconnell.copartisan <- -0.165048532 # copartisan (willing to take lower wages)
mcconnell.outpartisan <- 0.004853362 # outpartisan (demand higher wages)

#how much money are they willing to give up to work at a workplace democracy?
#find the ratio between copartisan and outpartisan
#copartisan
copartisan.ratio <- mcconnell.copartisan/copartisan.estimate
#outpartisan
outpartisan.ratio <- mcconnell.outpartisan/outpartisan.estimate

#multiply the ratio against our workplace democracy estimates to find willingness to pay
#copartisan estimate
esop.copartisan <- 0.2016158916*copartisan.ratio #ESOPs
codeterm.copartisan <- 0.1486489704*copartisan.ratio # codetermination
election.copartisan <- 0.0809887104*copartisan.ratio # manager elections

#outpartisan estimate
esop.outpartisan <- 0.2016158916*outpartisan.ratio #ESOPs
codeterm.outpartisan <- 0.1486489704*outpartisan.ratio # codetermination
election.outpartisan <- 0.0809887104*outpartisan.ratio # manager elections

#multiply against their dollar estimates
hourly.rate.sd <- 1.39 * 60/15.1 #this reflect the average work time in their study multiplied by the standard deviation

hourly.rate.sd*esop.copartisan*40*52
hourly.rate.sd*codeterm.copartisan*40*52
hourly.rate.sd*election.copartisan*40*52

hourly.rate.sd*esop.outpartisan*40*52
hourly.rate.sd*codeterm.outpartisan*40*52
hourly.rate.sd*election.outpartisan*40*52

#### ANALYZE DATA (NO CONTROLS) - under 80k ####
df <- df %>% filter(cje_income %in% sort(unique(df$cje_income))[1:6])
#Do people prefer firms with these features?
main_work_td <- lm_robust_omitted(work_binary ~ 
                                    copartisan +
                                    outpartisan +
                                    as.factor(cje_corporate_gov) +
                                    as.factor(cje_corporate_resp) +
                                    as.factor(cje_firm_size) +
                                    as.factor(cje_gender_ownership) + 
                                    as.factor(cje_healthcare) + 
                                    as.factor(cje_hours) + 
                                    as.factor(cje_training) + 
                                    as.factor(cje_location) + 
                                    as.factor(cje_sick_leave) + 
                                    as.factor(cje_parental_leave) + 
                                    as.factor(cje_race_ownership) + 
                                    as.factor(cje_retirement) + 
                                    as.factor(cje_income) +
                                    as.factor(cje_type_of_work) +
                                    as.factor(cje_union) +
                                    as.factor(cje_work_from_home) +
                                    as.factor(cje_team_work) +
                                    as.factor(cje_work_culture)
                                  ,
                                  data = df,
                                  clusters = df$rid
) %>%
  filter(term != "(Intercept)") %>%
  mutate(term_clean = factor(gsub(x = term,pattern = "as.factor",replacement = "")),
         estimate = ifelse(is.na(estimate), 0, estimate),
         std.error = ifelse(is.na(std.error), 0, std.error),
         conf.low.90 = estimate - 1.64*std.error,
         conf.high.90 = estimate + 1.64*std.error) 
main_work_td

#our estimates
esop.estimate <- 0.304657477 # ESOPs
codeterm.estimate <- 0.188254519 # codetermination
election.estimate <- 0.186789953 # workers elect managers
copartisan.estimate <- 0.093534444 # copartisan
outpartisan.estimate <- -0.283689399 #outpartisan

#mcconnell et al estimates
mcconnell.copartisan <- -0.165048532 # copartisan (willing to take lower wages)
mcconnell.outpartisan <- 0.004853362 # outpartisan (demand higher wages)

#how much money are they willing to give up to work at a workplace democracy?
#find the ratio between copartisan and outpartisan
#copartisan
copartisan.ratio <- mcconnell.copartisan/copartisan.estimate
#outpartisan
outpartisan.ratio <- mcconnell.outpartisan/outpartisan.estimate

#multiply the ratio against our workplace democracy estimates to find willingness to pay
#copartisan estimate
esop.copartisan <- esop.estimate*copartisan.ratio #ESOPs
codeterm.copartisan <- codeterm.estimate*copartisan.ratio # codetermination
election.copartisan <- election.estimate*copartisan.ratio # manager elections

#outpartisan estimate
esop.outpartisan <- esop.estimate*outpartisan.ratio #ESOPs
codeterm.outpartisan <- codeterm.estimate*outpartisan.ratio # codetermination
election.outpartisan <- election.estimate*outpartisan.ratio # manager elections

#multiply against their dollar estimates
hourly.rate.sd <- 1.39 * 60/15.1 #this reflect the average work time in their study multiplied by the standard deviation

hourly.rate.sd*esop.copartisan*40*52
hourly.rate.sd*codeterm.copartisan*40*52
hourly.rate.sd*election.copartisan*40*52

hourly.rate.sd*esop.outpartisan*40*52
hourly.rate.sd*codeterm.outpartisan*40*52
hourly.rate.sd*election.outpartisan*40*52
